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DESCRIPTION 



METHOD OF COMPUTING FIR FILTER COEFFICIENT AND PROGRAM 

FOR COMPUTING SAME 

5 

CROSS REFERENCE TO RELATED APPLICATION 

This application is based upon and claims the benefit of priority from the prior 
U.S. provisional Patent Application No. 60/418,313, filed on October 15, 2002, the 
entire contents of which are incorporated herein by reference. 

10 

TECHNICAL FIELD 

The invention relates to a computational method for impulse response 
coefficients of a universal maximally flat FIR Filter, and to a computational program for 
the same. 

15 

BACKGROUND ART 

More than thirty years have passed since the publication of the seminal paper of 
O. Herrmann on the design of maximally flat FIR digital filters [1]. Of the body of the 
related research that has been published since Herrmann's article, the most significant 
20 contribution to the theory of lowpass maximally flat FIR filters has been Baher's results 
on the design of filters with simultaneous conditions on amplitude and group-delay 
response [2]. In recent years, it has been recognized that Baher's filters form a unifying 
class of maximally flat filters. This class encompasses various types of seemingly 
distinct systems [3]. Calling this class of filters the universal maximally flat digital 
25 filters, the authors simplified Baher's formula for the transfer function and showed that 
lowpass filters of even or odd lengths, with linear or nonlinear phase response 
characteristics, as well as fractional delay and half-band filters may be readily obtained 
by setting the design parameters in an appropriate manner. 

Notwithstanding universal maximally flat digital filters have an explicit formula 
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for their transfer functions and, thus, may be readily designed using a computer algebra 
system, the simplified formulas provided in [3] are still unwieldy. A direct use of the 
formulas involves computation of binomial coefficients with possibly non-integer 
indices that appear as terms of the summand in a three-fold sum. Binomial coefficients 
5 with non-integer indices cannot be computed by the simple method of Pascal's triangle. 
Another common problem with computations involving the binomial coefficients is due 
to their tendency to introduce very large integers to the intermediate steps of the 
computation. With the flexible computer power today, this may not be a serious 
impediment in many situations. However, there can be cases where limitations on 
10 hardware and software resources dictate the use of efficient means in computation of 
filter coefficients. It is desirable that such methods be free of binomial coefficients and 
have a simple iterative structure. Moreover, it is known that if the design parameter 
associated with the group delay is limited to rational numbers, as is the case with many 
realistic settings, the values of the impulse response coefficients are rational numbers. In 
15 such cases, it is important to ensure that the answer to the numerical problem in 
computation of an impulse response coefficient is exactly a fraction, say 1/3, not a 
floating number that may be shown as "0.333333335". Therefore, it is desirable for a 
computation algorithm to maintain this property of the impulse response coefficients 
and yield rational values for rational design parameters. 
20 The following are prior art references related to the present invention. 

[1] O. Herrmann, "On the approximation problem in nonrecursive digital filter design," 
IEEE Trans. Circuit Theory, vol.CT-18, no.3, pp.411-413, May 1971. 
[2] H. Baher, "FIR digital filters with simultaneous conditions on amplitude and delay," 
Electron. Lett., vol.18, pp296-297, 1982. 
25 [3] S. Samadi, A. Nishihara and H. Iwakura, "Universal maximally flat lowpass FIR 
systems," IEEE Trans. Signal Processing, vol.48, pp. 1956-1964, 2000. 
[4] R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics. 
Addison-Wesley, 1989. 
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DISCLOSURE OF INVENTION 

An object of invention is to compute impulse response coefficients of the 
universal maximally flat FIR filter efficiently in a short period of time. 

According to one aspect of the present invention, the computer executes a first 
5 operation by a first recurrence formula, receiving a filter order (positive integer) of a 
universal maximally flat FIR filter, the number of zeros at z=-l (integer equal to or more 
than zero), and a parameter for a group delay at z=l (rational number). The recurrence 
formula includes parameters the filter order, the number of zeros, and the group delay, 
and provides coefficients in Bernstein form representation of a transfer function of a 
1 0 universal maximally flat FIR filter. 

The first recurrence formula, for example, is expressed as 

V = (-l){(2d) bj.r + ( j - 1 ) bj. 2 '} / (N - j + 1 ) where 1 ^j^N with b 0 ' = 1 
andb-i' = 0, 

in which the filter order is N, the parameter for the group delay is d, the coefficients in 
15 Bernstein form representation of a transfer function of the universal maximally flat FIR 
filter are bj'. The resultant of the first operation is expressed as 
B'={l,bi\...,b N -K\0,...,0} where the number of zeros is K. 

Using a resultant of the first operation as an initial value, the computer executes 
the second operation by the second recurrence formula composed of additions, 
20 subtractions, and divisions by 2 to extract impulse response coefficients of the 
universal maximally flat FIR filter. 

The second recurrence formula, for example, is expressed as 

hi (p) = ( 1 + E ) hi*-" / 2 + (1 - E) hi-i*" l) / 2 where 1 ^p^N, O^i^p with h 0 (0) 

= B' and h.i (p) ={0,...,0}, 
25 in which a sequence for computing impulse response coefficients of the universal 
maximally flat FIR filter is expressed as hi^Kkj^Mko^hu^V •-)» and an arbitrary 
sequence Ai is expressed as E j = E ( E^ 1 Ai ) x B l A { = EAi = Am N E°Ai = Ai in which a 
forward shift operator satisfying the expression is E. The impulse response 
coefficients extracted from the resultant of the second operation are expressed as 
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h i =h if o (N) where O^i^N. 

BRIEF DESCRIPTION OF DRAWINGS 

Fig. 1 illustrates steps of introducing a computational algorithm for impulse 
response coefficients according to the present invention. 

Fig. 2 illustrates steps of proving the proposition 1. 
Fig. 3 illustrates steps of proving the corollary 1. 

Fig. 4 illustrates an overview of the computational algorithm for the impulse 
response coefficients according to the present invention. 

Fig. 5 is a flowchart showing the computational algorithm for the impulse 
response coefficients. 

Fig. 6 is a program list of the computational algorithm for the impulse response 
coefficients executed by computer. 

BEST MODE FOR CARRYING OUT THE INVENTION 

The embodiments of the present invention will be described with reference to 
the drawings in the following. 

First, a description will be made on mathematical underpinning of the 
computational algorithm for the impulse response coefficients according to the present 
invention. Fig. 1 illustrates steps of introducing the computational algorithm for the 
impulse response coefficients. 

It was shown in [3] that the transfer function of the universal maximally flat 
digital filters, a family of filters identical to those proposed by Baher [2], is given by 

H N>K> d(z)= H b, ( I^)Mi^-) N -i, (1) 

0<j<N-K 

where the coefficients bj are characterized by the three-term recurrence 

jbj + ldb^-G-N-^b^^O, j>l, (2) 

with the initial values 
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b 0 = 1 and b_, = 0. (3) 



By converting (1) to the power form representation 

H N ,K,d<z)= 51 H k s~ k (4) 

0<k<N 

the impulse response coefficients hk become [3] 

h.- z E E k „ 0 n. « 

0<5<N-K0<p<k0<i<j 

Obviously, (5) is a computationally expensive formula. To ameliorate this 
situation, devised is an algorithm that exploits the combination of (1) and (2) to compute 
hkfor k= 0,...,N. 

First, (1) is written in the traditional Bernstein form by introducing the sequence 
bj' defined by 

. 7§r. 0<j<N, 
bl d ^l O (6> 



0, otherwise. 



Then obtained is 



H N ,K,d(z) - 2_ b ) ( ; 1 t— 2~~ 3 ( ~~T~~ } 
0 <i<N-K V / 



)N-i. (7) 

'111 z z. 

0<i<N-K 



This is the Bernstein form of the transfer function in the traditional sense. 
Although (7) has an additional binomial coefficient compared to (1), we will see later 
that the introduction of bj' results in considerable simplification. 

A recurrence for bj' is obtained by substituting (6) into (2) and dividing both 
sides byj(?)=£0. 

Thus, 



b,' + 2d0b^-(j-N-2)0b|_ 2 = O, 1<j< 



N. (8) 



The above recurrence can be simplified by invoking the definition of the 
binomial coefficients [4] and a few simple algebraic steps. It follows that 
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with the initial values 

b£ = 1 and bi t =a (10) 

Thus, for given values of N and d the values of bj' can be computed recursively by (9). 

The next step is to convert (7) to the power form (4) and obtain the impulse 
response coefficients. A direct expansion of the summand in (7) using the binomial 
theorem will introduce additional sums and new binomial coefficients. Hence, this is not 
a desirable action to take. The other alternative is to use a conversion identity well 
known in the mathematics literature and in the field of computer-aided geometric design. 
The identity is based on taking successive differences of the Bernstein coefficients. We 
state this identity as a proposition and provide a proof for it. 
(Proposition 1) 

Allow the Bernstein form of the polynomial p(x) = 2 o^i^N Pi x l to be given by 

pW- £ btf^xMi-x)"^ (ii) 

0<i<N \ / 

Then 

\/S>h 0i i = 0 i _ ! N ! (12) 

where the operator A is the forward difference operator. 

The forward difference operator is defined by the general recurrence 

A J b l = A(A^ 1 b i ) f (13) 

together with 

A 1 b t = Ab t = b i+ i - b t> A 0 b t = b t . (14) 

For instance, Abo = bi - b 0 and A 2 b 0 = b 2 - 2bi + bo. In this notation, the 
operator acts on the indices, written as subscripts, of the members of a sequence. To 
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prove the proposition, we also need the forward shift operator E. The effect of E on the 
index i of a sequence like bi is specified by 

Pb i =E(Ei- 1 b i ). (I 5 ) 

together with 

E 1 bi = Eb t = b l+1 , E° W == bi. (16) 

E is related to A by 

A-E-1. (17) 

(Proof of proposition 1) 

Fig. 2 illustrates steps of proving the proposition 1. Starting with the Bernstein 
form of p(x) and expand the term (1 - x) N "* using the binomial theorem to obtain 



pM- E Z t-'J^fflfVW (18) 

0<i<N 0<j<N-l V / \ ' / 



By the symmetry of the binomial coefficients and the trinomial revision 
property [4], obtained is 

On the other hand 

b t = E l b 0 . (20) 

Thus, it can be written as 

pw- z z (-') N - i (•)( N N _Ti j ) %N " ,E,, ' o • m 

0<i<N0<j<N-i \' / \ V 

The ranges for the indices of the two-fold sum above may be interchanged as 

0<1<N 0<i<N-l 0<j<N 0<l<j 

with no effect on the value of the sum. Now, N - j is substituted for j in the sum, and 
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obtained is 



p W _ £ L(-.) N --<^'( N N _,)( N ^l' > !-i')- N - 1N - i 'E l b 0 . (23) 
0<N-j<N0<i<j \ V \ v 'V 



3<N-j<N 0<i<j 

The right side above simplifies to 



0<j<N 0<l<j V / V / 



5 Note that 



o<t<j v / 



(25) 



It then follows that 



(26) 



p(x)= £ ^(^(E-iybo 

that is the power form representation of p(x). From the uniqueness of the power from 
10 representation of a polynomial and the relation between the operators E and A it 
follows that pi = (?) A { b 0 . (the end of the proof) 

Now presented is a main and last result in this preparatory phase before 
proceeding to the derivation of the algorithm. 
(Corollary 1) 

15 The transfer function H N ,K,d( z ) is given by 

H N ,K,d(*) - (^4n + 4^ >N ^ (27) 

(Proof of corollary 1) 

Fig. 3 illustrates steps of proving the corollary 1. 
Let 

z- 1 = 1 - 2x. (28) 

20 

Then, it can be written that 
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H NlK> d(x)= £ W^)xi(l-x) N ->. (29) 
0<j<N-K \'/ 

By the above proposition obtained is 

H N ,K,d(x)= Y. ("Wo* 1 . (30) 

By substituting for x in terras of z l in the above relation and expanding the 
summand using the binomial theorem, obtained is 

By the trinomial revision identity obtained is 

The sum indices can also be modified to write 

HN,K.a(,)= £ h^lT^fii^^. (33) 

0<i<N\ 1 / Z l<j V V 

Using p = j - i as an index for the second sum, obtained is 

hnmh- Y. L (V) 1 !™- (34) 

0<i<N \ / 0£p<N-i \ K / 

This simplifies to 

H N ,K,a(z)= Z f^(^) l (1+|) N ^. 

0<i<N \ / 

Invoking the binomial theorem further makes it possible to write 



(35) 



H N> K,d(*) = 0 + f ~ | z" 1 ) N b£. (36) 

By using the relation between A and E, (27) follows, (the end of the proof) 

Next a description will be made on details of the computational algorithm for 

the impulse response coefficients according to the present invention. 

The algorithm consists of two major stages. In the first stage, by running the 

recurrence (9) from j = luptoj=N-K, the following sequence is obtained. 
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10 

B' = (1,b 1 '...,b^_ K ,0,...,0). (37) 

The entries bj' of B' for N - K < j ^ N are set to zero. The next step is to 
expand the transfer function using (27). This can be done by noting that 

for arbitrary integer p > 0. Evaluation of both sides above in the power form yields 

^ K (p) z -i = (1+1 + l^V') £_*i V -"z-\ (39) 

where hi (p) is a sequence of coefficients defined by 

^ = (hgJ) = (h|g,hg ) ,...), (40) 
and operator E effects a forward shift on the index j. Thus, it is obtained that 

It follows that 



hM = lllh^ 1 ' + ^phfcA (42) 
The ranges for i and p are specified by 

1<p<N, 0<i<p. (43) 

15 The initial values are given by 

h< 0) =B', h^-CO.O,...)- (44) 

This completely specifies the recurrence we use in the second stage of our 
algorithm. After completing the Nth step, at the point when p has reached N, the value 
of the i-th impulse response coefficient is given by 



Hi = first term of sequence h[ N} , i = 0,...,N. (45) 



20 
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Fig. 4 illustrates an overview of the computational algorithm for the impulse 
response coefficients according to the present invention. In Fig. 4, each node represents 
a sequence hi (p) . The nodes marked by 0's represent zero sequences. The recurrence is 
illustrated in Fig.4 for N = 3 and an unspecified value of K. It can be visually verified 
5 that the recurrence has a triangular structure. 

Fig. 5 is a flowchart showing the computational algorithm for the impulse 
response coefficients according to the present invention. 

Receiving a filter order N (positive integer) of the universal maximally flat FIR 
filter, a number of zeros K (integer equal to or more than zero) at z=-l, and a parameter 
10 durational number) for a group delay at z=l, the computer computes impulse response 
coefficients according to the following steps. 

In step S10, all ofN+1 elements B'[j] (O^j^N) of a single-dimension array B' 
are initialized to zero. In step S20, all of N 3 elements r[p,i,j] (O^p^N, O^i^N, 0^j 
^N) of a three-dimension array r are initialized to zero. In step S30, an element B'[0] 
15 of the single-dimension array B' is set to 1. 

In step S40, every element of the single-dimension array B 1 is decided by 
changing in sequence an index j from 1 to N-K in a recurrence formula B 1 !]] = (-l)X 
{(2d)B , [j-l] + G-l)B f [j-2]} / (N - j + 1). That is, the operation by the recurrence formula 
(9) is executed. 

20 In step S50, elements r[0,0,j] (0^j^N-K) of the three-dimension array r are set 

to elements B'[j] (0 ^ j ^2 N-K) of the single-dimension array B\ 

In step S60, every element of the three- dimension array r is decided by 
sequentially changing, in the order of indexes j, i, p, an index j from zero to N-p, and an 
index i from zero to p, an index p from 1 to N in a recurrence formula r[p,i,j] = 
25 ( r[p-l,i-l,j] - r[p-l,i-l j+1] ) / 2 + ( r[p-l,i,j] + r[p-l,i,j+l] ) / 2. That is, the operation 
by the recurrence formula (42) is executed. 

In step S70, N+l elements h[i] (O^i^N) of a single-dimension array h are set to 
elements r[N,i,0] (O^i^N) of the three-dimension array r. 

In step S80, the single-dimension array h is output as impulse response 
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coefficients of the universal maximally flat FIR filter. This completes the 
computational processing of the impulse response coefficients by the computer. 

Fig. 6 illustrates a program list of the computational algorithm for the impulse 
response coefficients by computer. A pseudo-code of a procedure that returns the value 
5 of impulse response coefficients as an array is given in Fig.6. 

The following example helps elucidate the details. Let N = 3, K = 1 and d = 
(-l)/4. The actual values of the sequences hi (p) , i = 0,...,3, corresponding to the 3-D 
array r in the procedure of Fig.6, are as follows: 

h<* h" 

p=0: (1,1/6,-11/24,0) 

p = 1: (7/12,-7/48,-11/48) (5/12,5/16 -11/48) ( 46 ) 

v = 2 : (7/32, -3/16) (35/48, 1/12) (5/96, 13/48) 

p = 3: (1/64) (39/64) (31/64) (-7/64) 

10 Only non-zero values of the sequences hi (p) , i = 1,...,3, are presented. Note that 

all computations have been done in the field of rational numbers in an exact manner. 
Hence, round-off errors are avoided. 

As described above, an alternative representation for the transfer function of 
universal maximally flat FIR filters enables the development of a very simple algorithm 

15 for evaluation of the impulse response coefficients of the filters. The alternative form of 
the transfer function is based on the shift operator E and is derived from a technique 
referred to as finite calculus in [4]. 

The algorithm consists of two stages. In the first stage a two-term recurrence 
relation is used to obtain a 1-D sequence of numbers. If the value of the group delay 

20 parameter d is a rational number, then the outputs of the recurrence are rational. The 
outputs are then used in a simple N-step algorithm of triangular structure. In each step 
of this algorithm, the values of the sequences obtained in the preceding step undergo 
simple additions, subtractions and divisions by 2 to yield sequences that are shorter in 
size. Upon the completion of the Nth step, N sequences are obtained. The first entries of 
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these sequences are the values of the impulse response coefficients. As long as the 
computations in the second stage are done in the field of rational numbers, the final 
values remain rational This provides an accurate and fast means for exact computation 
of the impulse response coefficients. 

Here, a worst-case analysis of the computational complexity of the present 
algorithm will be described. 

The computational complexity may be analyzed by evaluating the complexity of 
the simple for-loop marked by (A) and the nested for-loop designated by (B) in Fig. 6. 
For the number of additions and subtractions C ± , it is obtained that 

c±= £ 4+ £ Z L 3 - (47) 

1 <j<N-K 1 <p<N 0<i<p 0<j<N-p 

For a closed-form expression, note that 

C±=4(N-K)+3 H (N-p + 1), (48) 

1<p<N 0<l<p+1 

that can be evaluated as 

C±=4(N-IC) + 3 £ (p + 2)(N-p). < 49) 

0<p<N 

Thus, obtained is 

C± = 4 (N - K) + 3 (| N + N 2 + 1 >4 3 ). (50) 

For the number of multiplications C* obtained is 

C*= 21 2 = 2(N-K). (51) 

1<J<N-K 

The number of divisions C/ is given by 

c /= Z 1+ Z 21 H 2 ' (52) 

1 <1<N-K 1 <p<N 0<i<p 0<j<N-p 

that becomes 

C, = (N - K) + 2(| N + N 2 + 1 N 3 ). (53) 
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In comparison, for the direct evaluation using (5), the number of additions and 
subtractions C ± Direct is . 

c Direct > ^ H H 1 (54) 

0<)<N-K 0<p<kO<i<> 

The right side above is of the order 0(N 4 ). The number of multiplications C* Durect , on the 
other hand, is 

cP irect > I I) k « ( 55 ) 

0<j<N-K 0<p<Tc 0<i<j 

The right side above is of the order 0(N 5 ). For the number of divisions C/ Direct in a direct 
evaluation, it can be written in 



-Direct 



> 21 H I> < 56 ) 

0<j<N-KO<p<kO<i<i 



10 The right side above, again, is of order OCN 4 ). Thus, in addition to the simplicity of 
most of the division operations in the present algorithm (that are divisions by 2), a 
significant reduction in the overall complexity is accommodated in terms of sheer 
numbers. 

The invention is not limited to the above embodiments and various 
15 modifications may be made without departing from the spirit and scope of the invention. 
Any improvement may be made in part or all of the components. 



INDUSTRIAL APPLICABILITY 

As described above, the impulse response coefficients of the universal 
20 maximally flat FIR filters may be computed efficiently using a two-stage algorithm. For 
a filter of order N with K zeros at z = -1, the first stage of the algorithm uses a two-term 
recurrence to generate a sequence of length N - K + 1. This sequence undergoes an 
N-step iterative process that involves additions, subtractions and divisions by 2. The 
algorithm does not involve any form of binomial coefficients. Moreover, the algorithm 
25 requires less number of arithmetical operations compared to a direct evaluation using 
the closed-form formula. The algorithm is mostly suitable for real-time generation of the 
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coefficients on a DSP chip because it may be easily implemented in a computational 
environment that has restrictions on the available hardware or software resources. 



